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The ability to quantify stochastic fluctuations present in biochemical and other systems is becoming increasing 
important. Analytical descriptions of these fluctuations are attractive, as stochastic simulations are compu- 
tationally expensive. Building on previous work, the linear noise approximation is developed for biochemical 
models with many compartments e.g. cells. The procedure is then implemented in the software package 
COPASI. This technique is illustrated with two simple examples and is then applied to a more realistic 
biochemical model. Expressions for the noise, given in the form of covariances matrices, are presented. 



PACS numbers: 87.18.Tt, 05.40.-a, 87.18.Vf 



I. INTRODUCTION 

When modeling chemical or biochemical reactions, 
most of the time the particulate nature of the con- 
stituents is ignored, and their concentrations are treated 
as deterministic and continuous variables. This is the 
standard thermodynamic approach to reaction kinetics. 
However there are a number of situations where this ap- 
proach is not validi"— . These situations usually arise in 
the case of biochemical reactions in a cellular environ- 
ment. Here, there may only be hundreds, or even fewer, 
of molecules of one particular type and the individual 
nature of the molecules cannot be disregarded. 

It might be thought that, since statistical fluctuations 
involving N particles typically have an effect of size 
, the usual thermodynamical description would 
be approximately valid, with errors of less than a few per 
cent for systems of this kind. This is not the case, as 
many simulations of biochemical reactions, and also ex- 
periments, have shown e.gj^i^. There are cases where a 
stationary state with constant concentrations is expected 
from an analysis of the deterministic rate equations, but 
simulations of the stochastic model show sustained cy- 
cles and other complex oscillatory behavior'''^. So, far 
from being only slightly quantitatively incorrect, the rate 
equations can give qualitatively wrong results. 

There are two principal methods of investigating the ef- 
fects of statistical fluctuations on reaction kinetics. The 
first is computational. When the reactions are formu- 
lated as a chemical master equation there are, formally 
exact, simulation algorithms which theoretically allow for 
a complete understanding of the nature of the resulting 
behavior— lis. However they are time consuming, and for 
many realistic biochemical reactions prohibitive in terms 
of computational time. Fortunately a second, analytic 
but approximate, approach is available which is usually 
accurate and easy to use. This is the linear noise ap- 
proximation (LNA) where stochastic deviations from the 
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deterministic description are kept to linear order only. 
The spectrum of these fluctuations can then be calcu- 
lated analytically^. 

Both of these approaches have been systematically for- 
mulated in the last few years. Stochastic simulations 
are one important feature of the freely available software 
COPASI^^'^^ which gives researchers interested in mod- 
eling and simulation of biochemical networks easy access 
to a variety of different analysis and simulation methods. 
The LNA may be formally derived from the chemical 
master equation by using the system-size expansion of 
van Kampeni, and several new insights have been ob- 
tained in this way. In particular, the reason why fluctu- 
ations are still important when TV is relatively large can 
be understood as an amplification phenomenon; a res- 
onance effect means that the size of the fluctuations is 
of the order of C/\/N, where the constant C is not of 
order one, but of order of 10^ or 10^ or even largeri^f^. 
Therefore fluctuations may be significant even when N 
is of the order of a million. 

In this paper we further develop the theory of the LNA, 
discuss it in the context of COPASI, and apply the results 
to specific examples of biochemical reaction networks. In 
a previous paper^^, we showed how to incorporate the 
LNA into the COPASI program. There we combined the 
LNA with optimization or parameter scanning in a closed 
loop to study the correlation structure of biochemical sys- 
tems and illustrated the use of this new strategy on sev- 
eral examples. However that study was restricted in an 
important way. In many situations molecules react in a 
region with a given volume (denoted as a compartment), 
but may also be able to move from one compartment to 
another, possibly undergoing a chemical transformation 
on the way. The volumes of these compartments may dif- 
fer, and these differences will have important effects on 
the effective concentration fluxes. In our previous work 
we avoided this complication by only investigating mod- 
els with a single compartment. Biochemical models can, 
in principle, be reformulated to refer to only one compart- 
ment. However this is a serious restriction, as this proce- 
dure must be done 'by hand' and is both time consuming 
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and prone to error. Thus, in this paper, we show how 
to generahze the analysis to multi-compartment models, 
and so considerably extend the number of systems which 
can be investigated using our methodology. 

The outline of the paper is as follows. In Sec. HI] we 
introduce the formalism of the LNA when two compart- 
ments are present in the model and illustrate its use on 
a simple example. We then generalize this treatment in 
Sec. mil to many-compartment systems and discuss how 
the method can be used within COPASI, using a slightly 
more complicated model to illustrate this. An applica- 
tion to a biochemical model taken from the literature is 
given in Sec. IIVI to show the power of the method. We 
conclude in Sec. V with a discussion of possible future 
avenues of investigation. The Appendices contain some 
mathematical details, as well as some specific details of 
the models analyzed in the main text. 



II. FORMALISM 



When discussing the LNA in detail, we shall use the 
same formalism as in our previous paper 1'^. The general 
biochemical model, with K species Yi , . . . 
reactions, will be written as 



Yf. and M 



riiYi 



Pki^k 



(1) 



^ku^k 



PimYi + 



PkM^i 



where the r^^ and pi^j {i = 1, . . . , K; fi — 1, . . . , M) are 
the numbers of molecules of reactants and products in- 
volved in each reaction. This may also be written as 



K 



K 



(2) 



The elements of the stoichiometric matrix, 
Tin , describe how many molecules of species Yi are gained 
or lost due to reaction fi. Biochemical systems frequently 
contain conservation relations. That is, although there 
are K chemical species present in the system, they may 
not all be able to vary independently. We denote the 
number of conservation relations by A and the number 
of independent species by K=K — A. As ini^, we work 
with a system of dimension K: information about the 
remaining species may be obtained from the conserva- 
tion equations. When describing the biochemical model 
stochastically, we begin by defining a chemical master 
equation, where the state of the system is described by 
n = the vector of particle numbers for 

species Yi , . . . , Yk ■ This equation describes the time evo- 
lution of P{n,t): the probability for the system to be 



found in the state n at time t. It has the formM 
^El^ ^ J2 T,in\n - v^)P{n - v^, t) 
'T^{n + u^\n)P{n,t), 



(3) 



where Tp(n'|n) is the transition rate from state n to n' 
due to reaction ^. The new state n' is defined by the 
stoichiometry of the reaction, — (i^i^, . . . , vk^)- To 
illustrate how the method applies to systems with more 
than one compartment we shall look at a simple system 
involving the diffusion of molecules between two com- 
partments. Molecules in the first compartment, of vol- 
ume Vi, are labeled D, molecules in the second compart- 
ment, volume 1^2, are labeled E. The system is summa- 
rized in Table m where ni is the number of D molecules, 
and n2 is the number of E molecules. There is an ex- 
tra complication to consider when defining the transition 
rates for a system with many compartments. That is, 
how do the transition rates of compartment-crossing re- 
actions scale with the volumes of the compartments? The 
rates of transitions from one compartment to another one 
should be dependent on properties of the contact surface 
between the compartments. These properties, such as 
the size of the diffusion area or the number of channels, 
will generally not scale with any compartment volume in 
practice. Rather, this relation will be different for differ- 
ent reaction systems, and will depend on the geometry 
of the overall system, and on how the compartments are 
connected. The transition rates for both reactions con- 
sidered here are proportional to the concentration of the 
substrate, and scale with the volume of the first com- 
partment. In this paper, we shall just consider transi- 
tion rates that scale linearly with volume in this way. 
We shall return to this issue in Sec. |Vl There is a con- 



TABLE I. A two-compartment reaction system 



A simple reaction system showing diffusion between two 
compartments. 





Reaction 


Stoichiometry 


Transition Rate 


1 

2 


D — > E 
E — > D 


^2 = (1,-1) 


Ti{n + vi\n) — kni 
T2{n + U2\n) = fc^nz 



served quantity present in the system: the total num- 
ber of molecules, u, is unaffected by each reaction i.e. 
ni -\- n2 = u. Therefore, the system has only one inde- 
pendent variable, so we can describe the system by either 
the number oi D or E molecules present. Here we choose 
to use D. We rewrite the transition rate for the second 
reaction as T2 = fc(Vi/V2)(u — rii), using the conserva- 
tion equation. The contributions from each reaction to 
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Eq. ([3]) may be written as 

Reaction 1: {E+^ - l)[kniP{n,t)], 
Reaction 2: (E^^ - l)[k--f{u - ni)P{n,t)], 

(4) 

where the E are step operators, e.g. E~^/(ni) — f{ni — 
1). To perform the van Kampen expansion^, we split 714 
into a deterministic part, the macroscopic concentration, 
X, and a stochastic part, ^: 

After the change of variables given in Eq. ([S|), the step 
operator may be written as 

The master equation is then rewritten in terms of the 
new variables and terms of the same order are equated. 
The first order terms recover the macroscopic descrip- 
tion of the system. That is, an ODE for the chemical 
concentration x 

— — k[a'^{u — x) — x], 

(7) 

where — V1/V2 and u = limvj^-^ocu/Vi. Thus, the 
volume ratio becomes an extra parameter of the prob- 
lem. The macroscopic equation has a stable fixed point^^, 
which we shall denote hy x = x* . The second order terms 
define a Fokker-Planck equation^ for n(^,i), the proba- 
bility distribution of the random variable ^, 

This equation is fully characterized by A and i?, which 
are functions of the macroscopic concentration of the 
chemical species, x{t). In this example it is found that 
they are 

B = kx + ka^{u — x). 

(9) 

The general form of A and B will be discussed in Sec. lIIII 
For models with K > 1, the Fokker-Planck equation 
will be multi-variate, and A and B will be square ma- 
trices of dimension K. A and B are often called the 
drift and diffusion matrices, respectively. When they are 
evaluated at the fixed point, x = x* , A and B become 
constant and 11 describes the fluctuations around the 
macroscopic steady state. From the Fokker-Planck equa- 



tion, information about the moments may be obtained. 
Of particular interest in this paper are the covariances, 
^ij — ((Ci ~ ^ (Cj)))- The matrix of covariances 

satisfies the following matrix equation^ 

AE + EA^ + B = 0, (10) 

which is known as the Lyapunov equation^^. For our 
one-dimensional example, E is just the variance of the 
random variable i.e. E = ((f — (^))(^ — {£,)))■ As dis- 
cussed in our previous work, we are more interested in the 
variance in terms of molecular numbers ni, rather than 
^. That is, the desired form of the variance, C, in this ex- 
ample is C = ((rii — (rii))(ni — (ni))). From our previous 
paper^- (ggg alsoii) it is possible to link the two quan- 
tities by C = ViE. The general form of the relationship 
between C and S in a model with many compartments 
will be discussed in Appendix |^ For one-dimensional 
problems, Eq. ([10]) can be trivially rearranged to solve 
for E, the variance of the fluctuations around the steady 
state. We can compare this result from theory with that 
obtained via simulation, using the Gillespie algorithni^, 
as implemented in COPASI. The value for the variance 
found by using the LNA was 360, in units of particle num- 
bers squared, compared to 360.0 (to one decimal place) 
obtained from simulation, averaged over 5000 time se- 
ries. The numerical values used for the parameters in 
this example are as follows: Vi = 0.3, V2 = 0.2 (both in 
picolitres), k = 0.2 and the total number of molecules in 
the system was 1500. 



III. GENERAL CASE 



As mentioned in the previous section, the A and B 
matrices are functions of the macroscopic concentrations. 
But the transition rates that define the master equation 
depend on the molecular populations, which are discrete 
quantities. So, to define the general form of A and B 
we must define a macroscopic quantity F^{x) that cor- 
responds to the transition rate Tp(n + i^^\n) for reac- 
tion fi. To do this, we simply make the replacement 
Tii/V'^^^ — > Xi, since \imY(i)^^{ni) /V'^'^^ becomes equal 
to Xi in the thermodynamic limit. Here, V^^^ denotes the 
volume of the compartment within which species i is lo- 
cated. In this paper all volumes written with bracketed 
superscripts will be defined in this way: volumes with 
subscripts define the actual volume of a compartment. 
The general forms of A and B look very similar to those 
given ini^, except for extra volume factors, which are 
picked up from the step operators, as shown in Eq. ([B]). 
In one compartment models these volume factors cancel. 
The general form for A and B for multi-compartment 
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models are 



A/ 



^ ^/\n^)vTI) 



(11) 



It is convenient to define the above matrices in terms 
of the F^{x), since these are quantities which COPASI 
aheady calculates. In COPASI they are called 'particle 
fluxes', and are calculated as part of the Steady State 
Task. 

In one compartment models the entries of the matrix 
A are found to be identical to the entries of the Jaco- 
bian of the macroscopic system evaluated at the fixed 
point. If the general form of the ODEs is taken to be 
dxi/dt = giix), then Aij — dgi/dxj. This is not the 
case in models with many compartments. However, the 
correct form of A for many compartment models can be 
found by applying a similarity transformation to the Ja- 
cobian, which from now on we will call A. The relation- 
ship is found to be 



SAS- 



S 



diag(VVW, VV^, . . . , VVW). 

(12) 

As COPASI is able to calculate the Jacobian, we will uti- 
lize this calculation, instead of performing an extra one. 
However, there is an additional complication here, since 
COPASI uses ODEs for the expectation of the number of 
molecules of the chemical species, rather than their con- 
centrations. The Jacobian calculated from the former, 
which we will denote A, is not identical to the one cal- 
culated from the latter, A. However, these two matrices 
are similar, and can be easily related to each other (see 
Appendix \X\ for details): 

A^S2AS^\, 52 = diag(F(i\t/(2\...,F(^)). (13) 

We can use this relation, along with the relation linking 
A with A, the desired matrix, to find the relation between 
A and A. It is i = SAS-\ 

We want to use these matrices to define a Lyapunov 
equation which can be solved to yield the covariances. 
One course of action would be to convert A to the desired 
form. A, and use this to solve Eq. (fTO)) for 5 and convert 
to C, the covariances in terms of particle numbers, using 
Eq. (|A5[) . It is slightly more straightforward, however, 
to define an equivalent Lyapunov equation, involving A 
and C, 



AC + CA^ + B = 0, 



(14) 



where, in order for the above equation to be equivalent to 
Eq. (fTU]) . we make the identification 13 = SBS. This is 
equivalent to defining B in Eq. (1111) without the square 



rooted volume factors in the denominator. This is the 
form displayed in COPASI. COPASI then solves this 
equation for C, using the Bartels-Stewart algorithn>i^. 
The equivalence of the two Lyapunov equations is proved 
in Appendix [XI 



TABLE II. The three-compartment reaction system 





Reaction 


Particle Flux 


1 


Ge — > Gi 


Fi{x) = kiXfiVs 


2 


2Gi 3G2 + 4G3 


F2{x) = k2xlVi 


3 


G2 — >; G5 


Fz{x) = ^3X2X5^2 


4 


Ga-^ 


Fi{x) = k^x^iV^ 


5 


G4 — > G5 


F^{x) = k^XiVi 


6 


G5 — > G4 


Fe{x) = kex^Vi 


7 


Gi — > Gq 


Fr{x) = kjxiVs 


8 


-^Gi 


Fsix) - ksVi 



We end this section with another example, this time 
with three compartments instead of two, to illustrate the 
natural extension of the method to an arbitrary number 
of compartments. This three compartment model has six 
species. Species Gi and G4 are located in compartment 
of volume Vi , G2 , G3 and G5 in compartment of volume 
V2 and Gq in compartment with volume V3. The num- 
ber of Gi molecules is denoted by ni , the number of G2 
molecules by ^2, and so on. The reactions are described 
in Table ini along with their particle fluxes. 

In reaction 3, the rate of degradation of species G2 is 
now modified by the concentration of G5 within the com- 
partment. We notice that species G4 and G5 are linearly 
dependent on each other, and we choose to eliminate G5. 
The matrices A and B are given in Appendix [Bl CO- 
PASI solves the Lyapunov equation for the 5 dimensional 
'reduced' system, then recovers the full 6 dimensional co- 
variance matrix using the conservation relations, see^^i^ 
for details. 

Table IIIII shows the numerical values obtained for the 
covariances. The parameter values were chosen to be 
ki = O.ls-i, k2 = 0.02pl#-is-i, fc3 = 0.1pl#-is-i, 
/c4 = 2s-\ /c5 = O.ls-i, /c6 = O.ls-i, fcr = O.ls"! 
and k% = 50#pl~^s~^, where # denotes particle num- 
bers. The compartment volumes are Vi = 8, V2 = 100 
and V3 = 72 (in picohtres). Results for species G5 are 
not given, but may be found from conservation consid- 
erations. The conservation relation for this system is 
UA + n^ — 1300 molecules. 



IV. A BIOCHEMICAL APPLICATION 

In this section we look at a model due to Kongas and 
van Beek22, that studies the role of creatine kinase in the 
heart by examining energy metabolism in cardiac muscle. 
It is a two-compartment system, with a cytoplasm and an 
intermembrane space. We shall label the volumes of these 
compartments as Vc and Vi respectively. In the original 
article the model is studied deterministically. Here, we 
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TABLE III. Covariances for the three-compartment system 

The covariances of the fluctuations around the steady 
state, in units of particle numbers squared. Results 
obtained from the LNA are compared with those from 
simulation (in brackets), via the Gillespie algorithm. 
10000 time series, each of length 10000 seconds with 
5000 samples, were generated. 





Gi 


G2 


G3 


(^4 


Gg 




229 (229) 


-48 (-48) 


-55 (-55) 


0(0) 


-17 (-17) 


G2 


-48 (-48) 


841 (841) 


592 (592) 


34 (34) 


-87 (-88) 


Gs 


-55 (-55) 


592 (592) 


846 (846) 


0(0) 


-68 (-68) 


G4 


0(0) 


34 (34) 


0(0) 


89 (89) 


0(0) 


Ge 


-17 (-17) 


-87 (-88) 


-68 (-68) 


0(0) 


2396 (2393) 



will reduce the model in size, without changing the vol- 
ume ratio used in the article, to study it stochastically. 
The model is described schematically in Figure 1 o^. It 
involves 5 chemical species, ADP, ATP, creatine (Gr), 
phosphocreatine (PCr) and inorganic phosphate (Pi). 
All of these metabolites are present in both compart- 
ments, so we have 10 variables. We use a subscript, i, to 
denote the species in compartment Vi . The reactions are 
as follows: 



ADPi + Pi, 




ATP, 




ATP, + Cn 




ADP, 


+ PCr, 


ATP + Cr 




PCr + 


ADP 


ATP 




ADP- 


f Pi 


Pl^ 




Pi 




Cn 




Cr 




ADP, 




ADP 




PCr, 




PCr 




ATP, 




ATP. 





(15) 

A SBML (Systems Biology Markup Language^!) im- 
plementation of the model is available from the BioMod- 
els Databas o^^'^^ . This file, which can be downloaded 
and then read by COPASl, corrects an error in the rate 
equations given in^*^. As for our previous models, not 
all of the 10 species can vary independently, as there are 
three conservation relations present. Hence, COPASl re- 
duces the dimensionality of the model to 7. The conser- 
vation relations, in terms of molecule numbers, are 

AT Pi + ADP, + ATP + ADP = Ci , 
Pi - ADP, + PCr + Pii - ADP + PCr, = C2, 
Cr + Cr, + PCr + PCr, ^ C3, 

(16) 

where Ci, G2 and G3 are integer constants. We chose the 
values Gi = 35000, G2 = 1204, C3 = 11743 and did not 



alter the reaction parameters given in^S. The values for 
Gi, G2 and G3 were chosen to speed up the numerical 
simulations of the system (which are extremely slow) by 
reducing the overall number of molecules in the system, 
whilst ensuring that each species did not get too close to 
the zero particle boundary. The deterministic model of 
the system has a unique steady-state, which is described 
in Table IIVI We calculated the covariances of the fluc- 
tuations around the steady state using the LNA Task in 
COPASl. Table |V] shows these results and compares the 
values with those obtained from numerical simulation. 
Only results for 7 of the 10 species present in the model 
are shown: the covariances for the other species may be 
obtained from the conservation equations. 

TABLE IV. The creatine kinase model 

Description of the biochemical system at the 
steady-state. Steady-state values are given in terms of 
particle numbers, and are rounded to the nearest 
integer. 



Compartment 


Species 


Steady-State Value 




ADP 


17490 




ATP 


10 




Cr 


5859 




PCr 


13 




Pi 


16501 




ADP, 


11819 




ATP, 


5681 


Vi 


Cr, 


5857 




PCr, 


15 




P^^ 


13984 



V. DISCUSSION 

In our previous worki^ we automated the LNA pro- 
cedure for one compartment models into the software 
package COPASl. This allowed the fluctuations around 
a steady-state to be found without employing time- 
consuming algebra. It also means that fluctuation anal- 
ysis may be calculated in a closed-loop with other tasks 
available in COPASl, e.g. optimization. In this paper we 
extended this framework to include multi-compartment 
models. The accuracy of the method was demonstrated 
with simple examples in Sec. |ll] and IIIII In Sec. IIV| a 
more realistic model was considered. 

Some issues remain which have not yet been resolved. 
Reaction kinetics for cross-compartment reactions need 
to be carefully considered. For convenience, in many 
published models cross-compartmental reactions are of- 
ten defined so that they scale with volume in the same 
way as reactions in the 'bulk' of the cell: this is the ap- 
proach we have taken in this paper. In reality this will 
not always be the mentioned in Sec. [TTl The 
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TABLE V. Covariances for the creatine kinase model 

Values for the covariances of the fluctuations around the 
steady state, in units of particle numbers squared. 
Results obtained from the LNA are compared with 
those found from numerical simulation (in brackets), via 
the Gillespie algorithm. 4000 time series, each of length 
1000 seconds with 10000 samples, were generated. The 
typical standard deviations associated with the 
simulation results are less than 1. The largest, for the 
variances of AD Pi and AT Pi, arc 1.8. 





ADP, 


ATP 


Cr, 


PCr 


P^^ 


ATP, 


Cr 


ADP, 


8407 
(8406) 


-4 
(-4) 


1 

(1) 




(0) 


1512 
(1511) 


-2252 
(-2252) 


-1 
(-1) 


ATP 


-4 
(-4) 


10 

(10) 




(0) 




(0) 


-5 
(-5) 


-2 
(-2) 




(0) 


Cn 


1 

(1) 




(0) 


2936 
(2936) 


-6 
(-6) 


8 

(7) 


-1 
(-1) 


-2922 
(-2922) 


PCr 




(0) 




(0) 


-6 
(-6) 


13 
(13) 


-7 
(-7) 




(0) 


-6 
(-6) 


Pl^ 


1512 
(1511) 


-5 
(-5) 


8 

(7) 


-7 
(-7) 


9467 
(9469) 


-2732 
(-2731) 


7 

(7) 


ATP, 


-2252 
(-2252) 


-2 
(-2) 


-1 
(-1) 




(0) 


-2732 
(-2731) 


4846 
(4846) 




(1) 


Cr 


-1 
(-1) 




(0) 


-2922 
(-2922) 


-6 
(-6) 


7 

(7) 




(1) 


2936 
(2936) 



LNA will work with any particular scaling that is univer- 
sal for the system. However, a consistent methodology 
for analyzing fluctuations in systems with a mixture of 
scalings present (e.g. some reactions scale with a cell's 
surface area and some with its volume) has not yet been 
proposed. Note that when describing the reaction sys- 
tems, no mention of geometry of the compartments has 
been made. Any geometrical considerations should be re- 
flected in the specific form of the transition rates chosen 
by the user. 

Care should be taken when using the LNA in certain 
circumstances. As we mentioned in^"^, the Gaussian as- 
sumption of the noise breaks down if the system is close to 
a boundary e.g. zero molecules. In addition to this, the 
LNA for multi-compartment models is technically only 
valid when the compartment volumes are of comparable 
order. This is because terms of order e.g. (l/yTT) and 
(l/-\/T^) are equated when performing the expansion of 
the master equation. 

Studying biochemical models in silico is an important 
method by which to study biochemical systems. As mod- 
els become more complicated (and therefore more realis- 
tic), more and more models will include many compart- 
ments. Therefore it is important to describe how the 
LNA can be used to study multi-compartment models. 
It is possible to transform multi-compartment models 
to a one compartment model, (as we did iiii^) but this 
must be done 'by hand' and is time consuming. With 
the framework described in this paper, models may be 



studied in their original form, using our extended imple- 
mentation of the LNA in COPASI. We believe that this 
will be a valuable tool in the study of stochastic effects 
in biochemical systems. 
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Appendix A: Mathematical details 

1. Covariance transformation 

Here we use a simple example to show how the two 
covariance matrices, C and S, are related. The number of 
molecules of species 1 is ni and the number of molecules 
of species 2 n2. Species 1 is located within a compartment 
of volume Vi , species 2 within a compartment of volume 
V2. The change of variables given in Eq. ([5]) leads to 

ni = Vixi + v^Ci, n2 = V2X2 + \/\^^2- (Al) 
The quantity C12 is defined as 

C12 = {{ni - {ni)){n2 - (^2))). 



(A2) 



From Eq. (|Al[) we can express the expectation values of 
ni and n2 to be 



(m) = Vixi, (712) = V2X2. 



(A3) 



Note that this is only true to the order up to which we 
are working, (see e.g.'*'^'* for details). Using the three 
equations above, we find that Cy = \/ViV2 , which 
is equal to \/Vi\^'E.ij , where S is the covariance matrix 
for the random variable ^. In general, , the covariance 
between species i and j in terms of particle numbers can 
be written as 



Ci- 



(A4) 



This relationship can also be expressed as a matrix trans- 
formation. For a general system, with K species, the 
relationship is 



C^SES, S ^ diag(VU«, VVi^, • . • , Vv(^)). (A5) 



2. Relation between A and A 

In Section IIIII we discussed the relationship between 
matrix A, calculated from the van Kampen expansion, 
and A, the form of the Jacobian calculated by COPASI, 
using ODEs for the expectation of the particle numbers. 
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Eq. ([13]) gives the relationship between A and A, the Ja- 
cobian calculated by using ODEs for the concentrations 
of the chemical species. We will highlight the differences 
between these quantities with a simple example, just con- 
sidering reactions 1 and 7 in Table [ill These reactions 
make the following contributions to the macroscopic rate 
equations, 



^jj- = ^{kiXQ - kjxi), 
at vi 



= {krxi - kixe). 

at 



From the above equations, the entries of A are found to 
be 



Aii^~/3^kj, Aie=f3^ki 
Aqi = kj, Aqq = -fci, 



(A7) 



where — V3/V1. We can rewrite Eq. (|A6p in terms of 
the (ni) and (tiq) 

= ki{ne} - kj—im), 



= K7-j-r(rii) - ki{ne). 



For the equations above, the elements of A are 

ill = -P'^kj, ii6 = ki, 
Aqi = p'^ki, Aqg = -fci. 



(A8) 



(A9) 



By considering how the Xi vary compared to the (n^) it 
is possible to find the following relation by inspection 



A.. — 4 . 

'J ~ I/O") 



(AlO) 



Again, this may be written as a matrix transformation. 
For a general system, 

A^S2AS^\ ^2 = diag(F(i\^(2),...,F(^)). (All) 

A similar relation may be found between A and A 



A^SAS-\ ^ = diag(VvW, 

(A12) 

Putting all this together, we find a relation between A 
and A 



A = SAS- 



(A13) 



(A6) 3. The equivalence of the Lyapunov equations 



Here we will prove the equivalence of Eq. (jlOl) and 
Eq. (jl4p . the two Lyapunov equations defined in the pa- 
per. We start by rewriting equation (jlOp as 

AS-'^SE + ESS-^A'^ + B = Q, (A14) 

where S has its usual form. Next we pre- and post- 
multiply by S: 

SAS-^SES + SESS-^A'^S + SBS ^0. (A15) 



Using Eq. jM]), Eq. (|Il3| and defining B = SBS we 
can reduce the above equation to 



AC + CA^ +B = 0, 



(A16) 



which recovers (|14p . the form of the Lyapunov equation 
solved by COPASI. 



Appendix B: Model details 

For the reaction system described in Table HIl the form 
of A and B are as below, where = V1/V2 and — 
V3/V1. The species are ordered Gi, G2, G3, G4 and Ge: 
species G5 having been eliminated due to conservation. 
The concentrations x are evaluated at the fixed point 
values X = X*. A and B are 



A = 



( -4k2Xi - krl3^ 
6ak2Xi 
8ak2Xi 


V krl3 





-fcaE 









-fc4 








a'^k3X2 







fci/3 \ 
' 



-fcl / 



(Bl) 



B = 



/ I3'^{kix6 + kjxi) + Ak2xl 
—Qk2ax\ 
—Sk2ax\ 


— kix^P — k-jXiP 



—6k2axi 
9k2a'^xl + k3X2T 
I2k2a^xl 





—%k2ax\ 
I2k2a'^xl 
l&k2a^x\ + knX'i 











k^X4^ -I- fcgF 




l3{kiXQ + hxi) \ 



(B2) 



kixe + k^xi J 
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where T = u — o?x\. 
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